function handles = currentFrequencyAnalysis(handles)
%
% traces: data array in which  rows are  samples and  cols are traces.
% Isteps: vector containing amplitudes of steps.
% startVal and endVal: indices (equivalently, sample number) for the 
%    beginning and end of the analysis window (i.e., the current step).
%    Defaults are startVal = 1000, endVal = 6000.
% threshVal: threshold for event detection, in mV. Default is threshVal = 5.
% sr: sample rate. Default is sr = 10000.
% 
% Typical usage (as of July 2005)
% - For a family of current step responses, first load the traces as 
%   trace type 'general physiology traces', load method 'selected traces'.
% - In Generic Trace Analysis, browse the traces and check which ones you 
%   would like to analyze (e.g. traces 1 to 11, or whatever).
% - Send the handles variable to the workspace, using the pushbutton.
% - Define the Isteps from the command line (see above)
% - Optional: define the startVal and endVal (in samples), threshVal (in mV), 
%   and sr -- if left out, these values will default to values listed above.
%
% Example:
% (First load the traces and send 'handles' variable to the workspace.)
% startVal = 1000; endVal = 6000; threshVal = 5; sr = 10000;
% traceVector = [1:7];
% traces = handles.data.map.mapActive.bsArray(:, traceVector); figure; plot(traces);
% Isteps = (0:100:1000);
% data = currentFrequencyAnalysis(traces, Isteps, startVal, endVal, threshVal, sr);
% 
% Display:
% First plot: traces
% Middle plot:
%
% Editing:
% adapted from Leopoldo's spikeFinder
% gs june 2005 -- added calculation
% gs july 2005 -- made it more user friendly and generic
% gs july 2005 -- incorporated into mapAnalysis
% ----------------------------------------------------------------------------
sr = handles.data.acq.sampleRate;
traces = handles.data.map.mapActive.bsArray(:, handles.data.analysis.currentFrequencyAnalysis.tracesToAnalyze);
Isteps = handles.data.analysis.currentFrequencyAnalysis.Isteps;
startVal = handles.data.analysis.currentStepStart;
startInd = round(startVal * sr);
currentStepDuration = handles.data.analysis.currentStepDuration;
endVal = startVal + currentStepDuration;
endInd = round(endVal * sr);
threshVal = handles.data.analysis.spikeThreshold;

[samples, numTraces] = size(traces);

if numTraces < size(Isteps,2)
    Isteps = Isteps(1:numTraces);
end

% stepDur = (endVal - startVal) / sr;

tracesOffset = [];
offsetVal = 0;
offsetAmount = 50;

for N = 1 : numTraces
    
    traceData = traces(startInd:endInd, N);
    filtered = medfiltLocalVersion(traceData, 100);
    traceFilt = traceData - filtered'; % gs20080308 had to transpose
%     traceFilt = traceData - medfilt1(traceData, 100);

    a = traceFilt > threshVal; % thresholding
    
    spikeRight = find(diff(a) < 0);
    spikeLeft = find(diff(a) > 0);
    n = length(spikeRight);
    
    if n == 0 
        spikeIndex(N) = NaN;
        spikeAmplitude(N) = NaN;
        spikeTimes(N) = NaN;
        meanSpikeRate(N) = 0;
    else
        for i = 1 : n
            [Y I] = max(traceData( spikeLeft(i) : spikeRight(i) ));
            spikeIndex(i) = I + spikeLeft(i) - 1;
            spikeAmplitude(i) = traceData( spikeIndex(i) );
            
            spikeTimes = spikeIndex / sr;
            
            meanSpikeRate(N) = round(n / currentStepDuration);
        end
    end
    
    if n >= 2
        interSpikeInterval = ( spikeTimes(2) - spikeTimes(1) );
        instRate1and2(N) = round(1 / interSpikeInterval);
    else
        instRate1and2(N) = NaN;
    end
    
    % store results
    data.trace(N).spikeTimes = spikeTimes;
    data.trace(N).spikeAmplitude = spikeAmplitude;
    data.trace(N).meanSpikeRate = meanSpikeRate;
    data.trace(N).instRate1and2 = instRate1and2;
    
    % save an offset version of the trace
    offsetVal = offsetVal + offsetAmount;
    tracesOffset(:,N) = traces(:,N) + offsetVal;
    
end

data.meanSpikeRate = meanSpikeRate;
data.instRate1and2 = instRate1and2;

disp(' ');
disp(['data.meanSpikeRate = [' num2str(meanSpikeRate) '];']);
disp(['data.instRate1and2 = [' num2str(instRate1and2) '];']);
disp(' ');

% figure
hfig = figure('Color', 'w', 'Position', [79   228   767   443]);

hplot = subplot(1,3,1);
numberOfSamples = size(tracesOffset, 1);
X = 0:(numberOfSamples-1);
time = X/sr;
hax = plot(time, tracesOffset);
axis tight;
title('Traces');
xlabel('Time (sec)');
ylabel('Voltage (mV)');

hplot = subplot(1,3,2);
hax = plot(Isteps, meanSpikeRate, 'bo-');
title('Mean spike rate');
xlabel('Current step amplitude (pA)');
ylabel('Mean spike rate (Hz)');

hplot = subplot(1,3,3);
hax = plot(Isteps, instRate1and2, 'bo-');
title('Initial instantaneous rate');
xlabel('Current step amplitude (pA)');
ylabel('Rate (Hz)');


% ========= local version of medfilt1 ===============
function m = medfilt1(s, w)
	if nargin == 1,
		w = 5;
	end
	
	s = s(:)';
	w2 = floor(w/2);
	w = 2*w2 + 1;

	n = length(s);
	m = zeros(w,n+w-1);
	s0 = s(1); sl = s(n);

	for i=0:(w-1)
            % changed code to do circular filtering
 		m(i+1,:) = [s0*ones(1,i) s sl*ones(1,w-i-1)];
%        m(i+1,:) = [s(end-i+1:end) s s(1:w-i-1)];
	end
	m = median(m);
	m = m(w2+1:w2+n);
